clear all
% This file simulates changes in the policies to incentivice the purchas
% eof EV.
%%%%%%%%%%%%%%%%%%%%%%%%%%% This version considers constant total sales
%%%%%%%%%%%%%%%%%%%%%%%%%%% across simulations to address only substitution
%%%%%%%%%%%%%%%%%%%%%%%%%%% between inside options
%first version 07/15/2021
%last version 10/04/2021
load('results_final')
%%%%%%%%% This function computes the Bounds for the fixed cost"%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

theta1=theta1_updated;
theta2=theta2_updated;
ob=length(p);
xi=xi_updated;
mc_update=p_nt;

% Calculate expected purged cost
mc_0=mc_update./(1+tariffs);
delta_rec=x_d*theta1_updated(1:k_D);

num_sim=100;

% Create supply matrix characteristics switching production location


policies={'none','trade tariff','tax exemption','both'}; % policy can be 'tax exemption' 'trade tariff' 'none'
results={};

%%
sim_res_ms=zeros(length(p),length(policies));
sim_res_p=zeros(length(p),length(policies));
for jj=1:length(policies)
    
    policy=policies{jj};
    new_sales_tax=ones(size(p))*0.19;
    new_tariffs=ones(size(p))*0.35;

    if strcmp(policy,'none')
        cf_salestax=sales_tax;
        cf_tariffs=tariffs;
    elseif strcmp(policy,'tax exemption')
        cf_salestax=new_sales_tax;
        cf_tariffs=tariffs;
    elseif strcmp(policy,'trade tariff')
        cf_salestax=sales_tax;
        cf_tariffs=new_tariffs;
    elseif strcmp(policy,'both')
        cf_salestax=new_sales_tax;
        cf_tariffs=new_tariffs;
    else
        print('enter valid policy')
    end

% set new price vector
    mc_bl2=mc_0.*(1+cf_tariffs);
    p_nt2=mc_bl2;
    p_cf=(p_nt2.*(1+cf_salestax))./1000;
    
    expmu1=expmu(theta2,v,x2,p_cf);
    sj_out=mktshares(delta,expmu1,QI,ind_market);
    
%     tic
%     [sj_out, p_out] = new_eq(p_cf,mc_bl2,theta2,ns,ind_market,x2,cf_salestax, delta, A1, QI,year);
%     toc
    sim_res_ms(:,jj)=sj_out;
    sim_res_p(:,jj)=p_cf;
end

%% Calcualte the expected baseline variable profit
%save('output_sim.mat','results','-v6')
%load('sim_res.mat')
%load 'output_sim.mat'
%
% none_res=results{1};
% tariff_res=results{2};
% tax_res=results{3};
% both_res=results{4};
% 
% sim_res_ms=[none_res(:,1),tariff_res(:,1),tax_res(:,1),both_res(:,1)];
% sim_res_p=[none_res(:,2),tariff_res(:,2),tax_res(:,2),both_res(:,2)];

%v_c=getranddraw(length(theta2)-1, ns);
%expmu1=expmu(theta2,v_c,x2,(p_nt.*(1+new_sales_tax))./1000);
%sj_ob=mktshares(delta,expmu1,QI,ind_market);

%sim_res_ms(year<2017,:)=repmat(sj_ob(year<2017),1,4);

%% Matching fiscal expenditure to no tariff excemption policy

tax_ava=find((sales_tax==0.05));
s_tar=0.35;
p_cf=(p_nt.*(1+sales_tax))./1000; % Baseline prices
mc_cf=mc_0; % Base line marginal cost
%s_mkt=sim_res_ms(:,1);

sum_ms=accumarray(ind_market,sj); % Total market share inside option
total_units=accumarray(ind_market,q); % Total units sold in each month
total_market=total_units./sum_ms;
sim_units=round(sim_res_ms.*total_market(ind_market),0);
[ww, zz] = ndgrid(ind_market,1:size(sim_units,2));

month_uni=accumarray([ww(:) zz(:)],sim_units(:));
month_ms=accumarray([ww(:) zz(:)],sim_res_ms(:));

total_market2=round(month_uni(:,1)./month_ms,0);

sim_units2=round(sim_res_ms.*total_market2(ind_market,:),0);

ST=sim_units2.*(sim_res_p-(sim_res_p./[1+sales_tax,1+sales_tax,1+new_sales_tax,1+new_sales_tax]));
TT=sim_units2.*(mc_0./1000.*[tariffs,new_tariffs,tariffs,new_tariffs]);
FC=ST+TT;
%p_ini=sim_res_p(:,2)./(1+sales_tax);
ss=2;
%%
options = optimset('Algorithm','sqp','GradObj', ...
'off','GradCon','off','DerivativeCheck','off', 'maxiter', 4000,...
'MaxFunEvals', 1000000,'TolFun',1e-12,'TolX',1e-12);

lb=0;%lower bound for mean and random coefficient for price = -inf
ub=0.35; %upper bound for mean and random coefficient for price = inf
% 
% for i=min(year):max(year)
%     tmp=year(year==1);
%     mc_cf_i=mc_cf(tmp,:);
%     ind_market_i=ind_market(tmp,:);
%     x2_i=x2(tmp,:);
%     [~,tax_ava_i]=ismember(tax_ava,tmp);
%     delta_i=delta(tmp,:);
%     QI_i=QI(tmp,:);
%     FC_C=FC(tmp,:);
%     
% 
s_tar2=0.19;


tic
[out_1,fval]=fmincon(@(ini2)new_eq_tt(ini2,mc_cf,theta2,ind_market,x2,tax_ava, delta, QI,FC,total_market2,ss,v,year),s_tar2,[], [], [], [], ...
lb,ub,[], options);
toc

cf_salestax2=ones(size(delta))*0.19;
tariff_cf2=(ones(size(delta))*0.35);
tariff_cf2(tax_ava)=out_1;
m0=mc_cf.*(1+tariff_cf2); 
p_out2=(m0.*(1+cf_salestax2));
expmu0=expmu(theta2,v,x2,p_out2./1000);
sj_out2=mktshares(delta,expmu0,QI,ind_market);

sim_un=round(sj_out2.*total_market2(ind_market,2),0);
st_0=sim_un.*(m0./1000.*(cf_salestax2));
tt_0=sim_un.*(mc_cf./1000.*(tariff_cf2));
total_0=tt_0+st_0;
f=abs(sum(FC(:,ss))*100-sum(total_0)*100);
%save('output.mat')
%% Carbon Tax
tier=zeros(size(p_nt));
UVT=9.37;

tier(p_nt<=1363*UVT,:)=0.015;
tier((1363*UVT<p_nt)&(p_nt<=3066*UVT),:)=0.025;
tier(p_nt>3066*UVT,:)=0.035;

tier(id_fuelcat=='Electric',:)=0.01;
tier(id_fuelcat=='Hybrid',:)=0.01;

qd=[0.02:0.02:1];

quan=quantile(emissions,qd);

[temp,cat] = histc(emissions,quan);
cat(cat==50,:)=49;
cat=(cat+1).*1.5* UVT;

m0=mc_cf.*(1+tariffs); 
p_out3=(m0.*(1+sales_tax+tier))+cat./1000;
expmu0=expmu(theta2,v,x2,p_out3./1000);
sj_out3=mktshares(delta,expmu0,QI,ind_market);



% tic
% [sj_out3, p_out3] = new_eq2(p_nt,mc_update,theta2,ns,ind_market,x2,sales_tax, delta, A1, QI,tier,cat,year);
% toc

% sj_out3(year<2019,:)=sj_ob(year<2019);
%save results_sim2
%% Carbon Tax 2
p_ini2=sim_res_p(:,1)./(1+sales_tax);

tier2=zeros(size(p_nt));
UVT=9.37;
ctt=1.5; % Multiplication parameter for emissions categories.

sc=2; % Scenario for comparation

qd=[0.02:0.02:1];

quan=quantile(emissions,qd);

[temp,cat2] = histc(emissions,quan);
cat2(cat2==50,:)=49;
cat2=(cat2+1);

lb=-inf;%lower bound for mean and random coefficient for price = -inf
ub=10; %upper bound for mean and random coefficient for price = inf

ini2=1;
outt=[0;0];
for kk=[2019,2020]
    tic
    [out_2,fval]=fmincon(@(tax) new_eq_tt2(tax,mc_cf,theta2,ind_market,x2, delta,...
      QI,FC,total_market2,ss,v,cat2,UVT,tariffs,sales_tax,tier,kk,year),ini2,[], [], [], [], ...
    lb,ub,[], options);
    outt(kk-2018)=out_2;
    toc
end

ctt_vec=ones(size(p)).*ctt;
ctt_vec(year==2019)=outt(1);
ctt_vec(year==2020)=outt(2);
cat_t=cat2.*ctt_vec* UVT;

m0=mc_cf.*(1+tariffs); 
p_out4=(m0./1000.*(1+sales_tax+tier))+(cat_t./1000);
expmu4=expmu(theta2,v,x2,p_out4);
sj_out4=mktshares(delta,expmu4,QI,ind_market);

sim_un=round(sj_out4.*total_market2(ind_market,2),0);
st_0=sim_un.*((m0./1000.*(sales_tax+tier))+(cat_t./1000));
tt_0=sim_un.*(mc_cf./1000.*(tariffs));
total_0=tt_0+st_0;

f=abs(sum(FC(year==2020,ss))*1000-sum(total_0(year==2020))*1000);
%save('CLCT.mat','sim_res_p','sales_tax','p_nt','emissions','mc_update','theta2','ns','ind_market','x2', 'delta', 'A1', 'QI','year','FC','total_market2','tariffs')

% tic
% [sj_out4, p_out4,ctt_c] = new_eq_ct(p_ini2,mc_update,theta2,ns,ind_market,x2,sales_tax, delta, A1, QI,year,FC,total_market2,tier2,cat2,UVT,sc,tariffs);
% toc

%ctt_c(year<2017,:)=1.5;
%sj_out4(year<2017,:)=sj_ob(year<2017);


%% plots variation in MS


% mean_sun=mean(tax_s(tax_s<0));
% sub=sales_tax;
% sub(sub==0.05)=mean_sun;
% p2=p_out2.*(1+sub);
% Get mean ms for fuel type for each scenario
sim_res_ms=[sim_res_ms,sj_out2,sj_out3,sj_out4];
sim_res_p=[sim_res_p, p_out2./1000, p_out3./1000,p_out4];

%sim_res_ms(year<2017,:)=repmat(sj_ob(year<2017),1,7);

tot_ms_sim=[];
fuel_cat=[];
s0_2=[];
mrkt_2=[];
mean_pr=[];

for qq=1:length(policies)+3
    temp=[];
    temp2=[];
    if qq==1
        for jj=1:max(ind_market)
            ttpp=accumarray(grp2idx(id_fuelcat(ind_market==jj)),sim_res_ms(ind_market==jj,qq));
            tempp=accumarray(grp2idx(id_fuelcat(ind_market==jj)),sim_res_p(ind_market==jj,qq),[],@mean);
            tot_ms_sim=[tot_ms_sim;ttpp]; 
            fuel_cat=[fuel_cat;unique(id_fuelcat(ind_market==jj))];
            s0_2=[s0_2;repmat(1-unique(s0(ind_market==jj)),size(ttpp))];
            mrkt_2=[mrkt_2;repmat(jj,size(ttpp))];
            mean_pr=[mean_pr;tempp];
        end
    else
        for jj=1:max(ind_market)
            ttpp=accumarray(grp2idx(id_fuelcat(ind_market==jj)),sim_res_ms(ind_market==jj,qq));
            tempp=accumarray(grp2idx(id_fuelcat(ind_market==jj)),sim_res_p(ind_market==jj,qq),[],@mean);
            temp=[temp;ttpp];
            temp2=[temp2;tempp];
        end
        tot_ms_sim(:,qq)=temp;
        mean_pr(:,qq)=temp2;
    end
end


day=ones(size(year));
date_mrk=unique(datetime(year,month,day));

plot(date_mrk,tot_ms_sim([1:4:288],1)./s0_2([1:4:288]))

hold on
plot(date_mrk,tot_ms_sim([1:4:288],2)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([1:4:288],3)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([1:4:288],4)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([1:4:288],5)./s0_2([1:4:288]))
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','fix fiscal cost','Location', 'southoutside')
title('Market Share Variation in Electric Cars for Different Policies')

plot(date_mrk,tot_ms_sim([4:4:288],1)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([4:4:288],2)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([4:4:288],3)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([4:4:288],4)./s0_2([1:4:288]))
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','Location', 'southoutside')
title('Market Share Variation in Internal Combustion Cars for Different Policies')

plot(date_mrk,tot_ms_sim([3:4:288],1)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([3:4:288],2)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([3:4:288],3)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([3:4:288],4)./s0_2([1:4:288]))
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','Location', 'southoutside')
title('Market Share Variation in Hybrid Cars for Different Policies')

 %% Plot variation price
plot(date_mrk,mean_pr([1:4:288],1))
hold on
plot(date_mrk,mean_pr([1:4:288],2))
hold on
plot(date_mrk,mean_pr([1:4:288],3))
hold on
plot(date_mrk,mean_pr([1:4:288],4))
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','Location', 'southoutside')
title('Price Variation in Electric Cars for Different Policies (1000 USD)')

plot(date_mrk,mean_pr([4:4:288],1))
hold on
plot(date_mrk,mean_pr([4:4:288],2))
hold on
plot(date_mrk,mean_pr([4:4:288],3))
hold on
plot(date_mrk,mean_pr([4:4:288],4))
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','Location', 'southoutside')
title('Price Variation in Internal Combustion Cars for Different Policies (1000 USD)')

plot(date_mrk,mean_pr([3:4:288],1))
hold on
plot(date_mrk,mean_pr([3:4:288],2))
hold on
plot(date_mrk,mean_pr([3:4:288],3))
hold on
plot(date_mrk,mean_pr([3:4:288],4))
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','Location', 'southoutside')
title('Price Variation in Hybrid Cars for Different Policies (1000 USD)')

%% Consumer surplus

CS_none = CS(theta2,v,x2,sim_res_p(:,1),delta,ind_market);
CS_tariff = CS(theta2,v,x2,sim_res_p(:,2),delta,ind_market);
CS_tax = CS(theta2,v,x2,sim_res_p(:,3),delta,ind_market);
CS_both = CS(theta2,v,x2,sim_res_p(:,4),delta,ind_market);
CS_fcf = CS(theta2,v,x2,sim_res_p(:,5),delta,ind_market);
CS_ct = CS(theta2,v,x2,sim_res_p(:,6),delta,ind_market);
CS_ct2 = CS(theta2,v,x2,sim_res_p(:,7),delta,ind_market);

CS_table=[CS_none, CS_tariff, CS_tax,CS_both,CS_fcf,CS_ct,CS_ct2].*1000;

plot(CS_none)
hold on
plot(CS_tariff)
hold on
plot(CS_tax)
hold on
plot(CS_both)
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','Location', 'southoutside')
title("Consumers' Surplus Variation")

save results_sim5
%% 

% plot(date_mrk,tar_s)
% hold on
% plot(xlim,[nanmean(tar_s(date_mrk>='01-Jan-2019')) nanmean(tar_s(date_mrk>='01-Jan-2019'))], 'r')
% legend('sale tax','mean sale tax','Location', 'southoutside')
% title("Counterfactual sale tax to match observed market share")
% hyp_tar=tariffs;
% hyp_tar(hyp_tar<0.35)=nanmean(tar_s(date_mrk>='01-Jan-2019'));

% p_cf2=(p_out2.*(1+new_sales_tax));
% mc_bl3=mc_0.*(1+hyp_tar);
%     
% tic
% [sj_out22, p_out22] = new_eq(p_cf2,mc_bl3,theta2,ns,ind_market,x2,new_sales_tax, delta, A1, QI,year);
% toc
% 
% sim_res_ms(:,5)=sj_out22;
% sim_res_p(:,5)=p_out22;

%% Tables

sim_units=round(sim_res_ms.*total_market(ind_market),0);
[ww, zz] = ndgrid(ind_market,1:size(sim_units,2));
[xx, yy] = ndgrid(year-2014,1:size(sim_units,2));

year_ms=accumarray([xx(:) yy(:)],sim_res_ms(:));
year_uni=accumarray([xx(:) yy(:)],sim_units(:));
month_uni=accumarray([ww(:) zz(:)],sim_units(:));
month_ms=accumarray([ww(:) zz(:)],sim_res_ms(:));

total_market2=round(month_uni(:,1)./month_ms,0);

sim_units2=round(sim_res_ms.*total_market2(ind_market,:),0);



%tot_mrk_tab=kron(total_market,ones(4,7)); % Market size table (4 scenarios, 4 fueltypes)
% tot_num_sim=round(tot_ms_sim.*tot_market2,0); % Table number of vehicles by scenario and by fueltype
% 
year_table=kron(unique(year),ones(48,1));
fuel_type_tab=kron(ones(72,1),[1:4]');

num_year_fc=zeros(4,7,2); %table with num of vehicles per year
for qq=1:7
    Num_unit_year=accumarray([year-2014 grp2idx(id_fuelcat)],sim_units2(:,qq));
    num_year_fc(:,qq,1)=Num_unit_year(end-1,:)';
    num_year_fc(:,qq,2)=Num_unit_year(end,:)';
end


Emi_unit_year=accumarray([year-2014 grp2idx(id_fuelcat)],emissions,[],@mean); % Average emissions by year and fuel cat


price_year_fc=zeros(4,4,2); %table with num of vehicles per year
for qq=1:6
    P_unit_year=accumarray([year_table fuel_type_tab],mean_pr(:,qq),[],@mean);
    price_year_fc(:,qq,1)=P_unit_year(end-1,:)';
    price_year_fc(:,qq,2)=P_unit_year(end,:)';
end

% Yearly market share sales

scenarios_ms=num_year_fc./sum(num_year_fc);% ms per scenario

% Data for Tables 12
Market_share_electric_row=[scenarios_ms(1,:,1);scenarios_ms(1,:,2)];
Market_share_hybrid_row=[scenarios_ms(3,:,1);scenarios_ms(3,:,2)];
 
sales=q.*p;

yearly_q=accumarray([year-2014,grp2idx(id_fuelcat)],sales);
yearly_tq=repmat(accumarray(year-2014,sales),1,4);

yearly_ms=yearly_q./yearly_tq;

%subplot(2,1)
bar(unique(year),yearly_ms(:,[1,3]))
legend('Electric','Hybrid', 'Location', 'southoutside')
xlabel('Year')
ylabel('Market sahre')

% Yearly q sold
yearly_q=accumarray([year-2014,grp2idx(id_fuelcat)],q);
yearly_tq=repmat(accumarray(year-2014,q),1,4);

yearly_q=yearly_q./yearly_tq;

subplot(2,1,1)
char1=bar(unique(year),yearly_ms(:,[1,3]))
xlabel('Year')
ylabel('Market sahre')

subplot(2,1,2); 
char2=bar(unique(year),yearly_q(:,[1,3]))
xlabel('Year')
ylabel('Percentage of units sold')


hL = legend([char1],{'Electric','Hybrid'});
% Programatically move the Legend
newPosition = [0.9 0.49 0.1 0.05];
newUnits = 'normalized';
set(hL,'Position', newPosition,'Units', newUnits);


% Mean price
yearly_p=accumarray([year-2014,grp2idx(id_fuelcat)],p,[],@mean);
yearly_p(yearly_p==0)=mean(yearly_p([1,3],4));

plot(unique(year),yearly_p(:,[1,3,4]))
xticks(unique(year))
legend('Electric','Hybrid','Gasoline', 'Location', 'southoutside')
xlabel('Year')
ylabel('MSRP + taxes (1000 USD)')
%% Fiscal cost computation
hyp_tar=new_tariffs; % hypothetic tariff for fixed fiscal cost counterfactual
sel_mrkt=ind_market(year>=2017);
for ii=min(sel_mrkt):max(sel_mrkt)
    temp=find(ind_market==ii);
    tar_1=hyp_tar(temp);
    [~,tax_ava_i]=ismember(tax_ava,temp);
    tax_ava_i(tax_ava_i==0)=[];
    tar_1(tax_ava_i)=out_1;
    hyp_tar(temp)=tar_1;
end

%hyp_tar(hyp_tar<0.35)=nanmean(tar_s(date_mrk>='01-Jan-2019'));



sim_res_pn=(mc_0./1000.*[1+tariffs,1+new_tariffs,1+tariffs,1+new_tariffs,1+hyp_tar,1+tariffs,1+tariffs]);
tariff_mat=(mc_0./1000.*[tariffs,new_tariffs,tariffs,new_tariffs,hyp_tar,tariffs,tariffs]);
sale_tax_mat=(sim_res_pn.*[sales_tax,sales_tax,new_sales_tax,new_sales_tax,new_sales_tax,sales_tax+tier,sales_tax+tier]);

sim_units=round(sim_res_ms.*total_market(ind_market,:),0);
sim_units(:,2)=round(sim_res_ms(:,2).*total_market2(ind_market,2));
sim_units(:,5)=round(sim_res_ms(:,5).*total_market2(ind_market,2));
sim_units(:,7)=round(sim_res_ms(:,7).*total_market2(ind_market,2));

% sale_tax_mat=(sim_res_p-(sim_res_p./[1+sales_tax,1+sales_tax,1+new_sales_tax,1+new_sales_tax,1+new_sales_tax,1+sales_tax,1+sales_tax]));
% sale_tax_mat(:,7)=(((sim_res_p(:,7)-(cat_t2./1000))./(1+sales_tax)).*(sales_tax));
% sale_tax_mat(:,6)=(((sim_res_p(:,6)-(cat./1000))./(1+sales_tax+tier)).*(sales_tax+tier));
variation_saletax=sale_tax_mat-sale_tax_mat(:,1);

[xx, yy] = ndgrid(year-2014,1:size(variation_saletax,2));
income_saletax=accumarray([xx(:) yy(:)],variation_saletax(:));
%tariff_mat=(mc_0.*[tariffs,new_tariffs,tariffs,new_tariffs,hyp_tar,tariffs,tariffs]);

other_tax=zeros(size(tariff_mat));
other_tax(:,6)=cat./1000;% Tax as stated in Colombian legislation proposal


other_tax(:,7)=cat_t./1000;% Counterfactual carbon tax coefficient (assuming no addicional tax on prices)


variation_tariffs=tariff_mat-tariff_mat(:,1);
income_tariffs=accumarray([xx(:) yy(:)],variation_tariffs(:));



fiscal_cost=sim_units.*(tariff_mat+sale_tax_mat+other_tax)*1000;
fiscal_cost_tab=accumarray([xx(:) yy(:)],fiscal_cost(:));
fiscal_cost_dif=(fiscal_cost_tab-fiscal_cost_tab(:,4));

% Data for table 12
fiscal_cost_table=fiscal_cost_tab(5:end,:);


%% Change in emissions
ave_km=20000; % assuming a total driving of 20k km per year

%sim_units=round(sim_res_ms.*total_market2(ind_market,:),0);
tot_emi=(emissions.*sim_units2.*ave_km)/(10^6);
life_emi=tot_emi*10; % 10 years Lifetime emissions in metric tons;
variation_emi=life_emi-life_emi(:,4);
yearly_emi_var=accumarray([xx(:) yy(:)],variation_emi(:)); % Emissions in metric tons

yearly_emi=accumarray([xx(:) yy(:)],life_emi(:)); % Emissions in metric tons
yearly_emi_dif=yearly_emi-yearly_emi(:,4);



cost_emi=fiscal_cost_dif./yearly_emi_dif; % cost in USD por ton of CO2
cost_emi(isnan(cost_emi))=0;
var_cost_emi=cost_emi-cost_emi(:,1);

cost_emi_tab=fiscal_cost_tab./yearly_emi;
cost_emi_dif=cost_emi_tab-cost_emi_tab(:,4);


% Profit per vehicle
% p_notax=(sim_res_p./[1+sales_tax,1+sales_tax,1+new_sales_tax,1+new_sales_tax,1+new_sales_tax,1+sales_tax,1+sales_tax]);
% p_notax(:,7)=((sim_res_p(:,7)-(cat_t2./1000))./(1+sales_tax));
% p_notax(:,6)=((sim_res_p(:,6)-(cat./1000))./(1+sales_tax+tier));
% 
% profit_v=(p_notax-(mc_0.*[1+tariffs,1+new_tariffs,1+tariffs,1+new_tariffs,1+hyp_tar,1+tariffs,1+tariffs]))*1000;
% 
% sim_units2=round(sim_res_ms.*total_market(ind_market,:),0);
% profit_tot=sim_units.*profit_v;
% profit_tot2=sim_units2.*profit_v;
% profit_avg=accumarray([xx(:) yy(:)],profit_v(:),[],@mean); % Profit per vehicle
% 
% year_uni2=accumarray([xx(:) yy(:)],sim_units(:));
% var_profit_avg=profit_avg-profit_avg(:,4);
% PR_yearly=accumarray([xx(:) yy(:)],profit_tot(:)); % Total Profit per vehicle
% 
% 
% PR_yearly2=accumarray([xx(:) yy(:)],profit_tot2(:)); % Total Profit per vehicle
% var_PR_yearly=PR_yearly-PR_yearly(:,4);
% var_PR_yearly2=PR_yearly2-PR_yearly2(:,4);

% Data for table 12
profit_avg_table=zeros(size(fiscal_cost_table));
profit_tot_table=zeros(size(fiscal_cost_table));
%% Consumer surplus
tot_units=accumarray([xx(:) yy(:)],sim_units(:));

% Welfare
year_month=accumarray(ind_market,year,[],@mean);
[zz, vv] = ndgrid(year_month-2014,1:size(CS_table,2));
[gg, hh] = ndgrid(grp2idx(id_market),1:size(sim_units,2));

CS_table2=CS_table.*1000;
units_months=accumarray([gg(:) hh(:)],sim_units(:));

CS_yearly=accumarray([zz(:) vv(:)],CS_table(:),[],@mean); % yearly consumer surplus



var_CS_yearly=CS_yearly-CS_yearly(:,4);
Total_CS=(CS_yearly.*tot_units);
vsr_TS=Total_CS-Total_CS(:,4);
%CS_per_con=CS_yearly./tot_units;

% Data for table 12
CS_yearly_table=CS_yearly(5:end,:);
Total_CS_table=Total_CS(5:end,:);

total_wel=(CS_yearly.*tot_units); % yearly welfare
total_wel_2=total_wel+fiscal_cost_dif;
total_wel_var=total_wel-total_wel(:,4);
total_wel_var_2=total_wel_2-total_wel_2(:,4);
%var_total_wel=total_wel-total_wel(:,1); % Variation in total welfare

% Data for table 12
total_wel_table=total_wel_2(5:end,:);

wel_emi=total_wel_var_2./yearly_emi_dif;

% weigthed average emissions per scenario


temp_emi=sim_units.*emissions;


fuel_type_tab=kron(ones(72,1),[1:4]');

emi_year_fc=zeros(4,7,2); %table with num of vehicles per year
for qq=1:7
    Emi_unit_year=accumarray([year-2014 grp2idx(id_fuelcat)],temp_emi(:,qq));
    emi_year_fc(:,qq,1)=Emi_unit_year(end-1,:)';
    emi_year_fc(:,qq,2)=Emi_unit_year(end,:)';
end

num_year_fc2=zeros(4,7,2); %table with num of vehicles per year
for qq=1:7
    Num_unit_year2=accumarray([year-2014 grp2idx(id_fuelcat)],sim_units(:,qq));
    num_year_fc2(:,qq,1)=Num_unit_year2(end-1,:)';
    num_year_fc2(:,qq,2)=Num_unit_year2(end,:)';
end

emi_year_type=emi_year_fc./num_year_fc2; % Weigthed average of emissions by vehicle type

tot_marlk_e=repmat(sum(num_year_fc2),4,1,1);

ms_year_emissions=num_year_fc2./tot_marlk_e;



% Market share and emissions by price
qd2=[0.2:0.2:1];

quan2=quantile(p,qd2);

[temp,p_cat] = histc(p,quan2);

p_cat=p_cat+1;

pr_year_fc=zeros(6,7,2); 
for qq=1:7
    pr_unit_year=accumarray([year-2014 p_cat],temp_emi(:,qq));
    pr_year_fc(:,qq,1)=pr_unit_year(end-1,:)';
    pr_year_fc(:,qq,2)=pr_unit_year(end,:)';
end

num_year_fc3=zeros(6,7,2); 
for qq=1:7
    Num_unit_year3=accumarray([year-2014 p_cat],sim_units(:,qq));
    num_year_fc3(:,qq,1)=Num_unit_year3(end-1,:)';
    num_year_fc3(:,qq,2)=Num_unit_year3(end,:)';
end

emi_year_p=pr_year_fc(1:end-1,:,:)./num_year_fc3(1:end-1,:,:); % Weigthed average emissions by price quintile

tot_marlk_p=repmat(sum(num_year_fc3),6,1,1);

ms_year_price=num_year_fc3./tot_marlk_p;

num_year_sg=zeros(8,7,2); 
for qq=1:7
    Num_unit_year_sg=accumarray([year-2014 grp2idx(id_staseg)],sim_units(:,qq));
    num_year_sg(:,qq,1)=Num_unit_year_sg(end-1,:)';
    num_year_sg(:,qq,2)=Num_unit_year_sg(end,:)';
end

emi_year_sg=zeros(8,7,2); 
for qq=1:7
    Num_emi_year_sg=accumarray([year-2014 grp2idx(id_staseg)],temp_emi(:,qq));
    emi_year_sg(:,qq,1)=Num_emi_year_sg(end-1,:)';
    emi_year_sg(:,qq,2)=Num_emi_year_sg(end,:)';
end
emi_year_sg=emi_year_sg./num_year_sg;

tot_marlk_sg=repmat(sum(num_year_sg),8,1,1);

ms_year_sg=num_year_sg./tot_marlk_sg;

% Mean profit

% temp_pi=sim_units.*profit_v;
% 
% num_year_fc4=zeros(4,7,2); %table with num of vehicles per year
% for qq=1:7
%     Num_unit_year4=accumarray([year-2014 grp2idx(id_fuelcat)],temp_pi(:,qq));
%     num_year_fc4(:,qq,1)=Num_unit_year4(end-1,:)';
%     num_year_fc4(:,qq,2)=Num_unit_year4(end,:)';
% end
% 
% profit_year_type=num_year_fc4./num_year_fc2; % Weigthed average of emissions by vehicle type

%% Individual consumer surplus

CS_none_ind = CS_ind(theta2,v,x2,sim_res_p(:,1),delta,ind_market);
CS_tariff_ind = CS_ind(theta2,v,x2,sim_res_p(:,2),delta,ind_market);
CS_tax_ind = CS_ind(theta2,v,x2,sim_res_p(:,3),delta,ind_market);
CS_both_ind = CS_ind(theta2,v,x2,sim_res_p(:,4),delta,ind_market);
CS_fcf_ind = CS_ind(theta2,v,x2,sim_res_p(:,5),delta,ind_market);
CS_ct_ind = CS_ind(theta2,v,x2,sim_res_p(:,6),delta,ind_market);
CS_ct2_ind = CS_ind(theta2,v,x2,sim_res_p(:,7),delta,ind_market);

CS_table_ind=[CS_none_ind, CS_tariff_ind, CS_tax_ind,CS_both_ind,CS_fcf_ind,CS_ct_ind,CS_ct2_ind]*1000;

temp_cs=sim_units.*CS_table_ind;

cs_year=zeros(4,7,2); %table with num of vehicles per year
for qq=1:7
    Num_cs_year=accumarray([year-2014 grp2idx(id_fuelcat)],temp_cs(:,qq));
    cs_year(:,qq,1)=Num_cs_year(end-1,:)';
    cs_year(:,qq,2)=Num_cs_year(end,:)';
end

cs_year_type=cs_year./num_year_fc2; % Weigthed average of emissions by vehicle type

cs_year_fc=zeros(6,7,2); 
for qq=1:7
    cs_unit_year=accumarray([year-2014 p_cat],temp_cs(:,qq));
    cs_year_fc(:,qq,1)=cs_unit_year(end-1,:)';
    cs_year_fc(:,qq,2)=cs_unit_year(end,:)';
end

save('final_data')

%% Expenditure

%exp_mat=sim_res_p.*sim_res_ms.*1000.*total_market2(id_market);
exp_mat=sim_res_p.*1000.*sim_units;
exp_year=accumarray([xx(:) yy(:)],exp_mat(:));
exp_diff=exp_year-exp_year(:,4);

%% Tables for paper

table_12_Row_Names={'Electric_vehicle_market_share_percent','Hybrid_vehicle_market_share_percent','Average_consumer_surplus_USD_per_consumer','Average_profit_USD_per_vehicle','Total_consumer_surplus_million USD','Total_profit_million_USD_per_market', 'Fiscal_cost_million_USD','Total_welfare_million USD per market'};
table_12_Var_Names={'Sales_tax_reduction_diff','Tariff_reduction_const_fical_cost_diff','Proposed_carbon_tax_const_fiscal_cost_diff'};

% Table 12 year 2019
format shortG
Table_12_2019=[
Market_share_electric_row(1,:)*100;
Market_share_hybrid_row(1,:)*100;
CS_yearly_table(1,:)/1000;
profit_avg_table(1,:);
Total_CS_table(1,:)/1.0e+6;
profit_tot_table(1,:)/1.0e+6;
fiscal_cost_table(1,:)/1.0e+6
total_wel_table(1,:)/1.0e+6;
];
Table_12_2019_diff=Table_12_2019-Table_12_2019(:,4);

Table_12_2019_final=[Table_12_2019_diff(:,3),Table_12_2019_diff(:,[5,7])];

Table_twelve_2019=array2table(Table_12_2019_final);
Table_twelve_2019.Properties.VariableNames=table_12_Var_Names;
Table_twelve_2019.Properties.RowNames=table_12_Row_Names;
writetable(Table_twelve_2019,'Table_12_year_2019.txt','Delimiter',' ','WriteRowNames',true)  

% Table 12 year 2020
format shortG
Table_12_2020=[
Market_share_electric_row(2,:)*100;
Market_share_hybrid_row(2,:)*100;
CS_yearly_table(2,:)/1000;
profit_avg_table(2,:);
Total_CS_table(2,:)/1.0e+6;
profit_tot_table(2,:)/1.0e+6;
fiscal_cost_table(2,:)/1.0e+6
total_wel_table(2,:)/1.0e+6;
];
Table_12_2020_diff=Table_12_2020-Table_12_2020(:,4);

Table_12_2020_final=[Table_12_2020_diff(:,3),Table_12_2020_diff(:,[5,7])];

Table_twelve_2020=array2table(Table_12_2020_final);
Table_twelve_2020.Properties.VariableNames=table_12_Var_Names;
Table_twelve_2020.Properties.RowNames=table_12_Row_Names;
writetable(Table_twelve_2020,'Table_12_year_2020.txt','Delimiter',' ','WriteRowNames',true)  


